/*
    functions for gridding temperature data

	This doesn't use a grid containing squares defined by latitude and
	longitude, as seems to be the standard, but instead distributes
	grid cells evenly across the globe.  This avoids any bias due
	to unequal land surface areas and also the problem of
	singularities at the poles.

    Copyright (C) 2012-2013 Bob Mottram
    bob@sluggish.dyndns.org

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "globalgrid.h"

globalgrid::globalgrid() {
}

globalgrid::~globalgrid() {
}

void globalgrid::threed_to_latlong(float x, float y, float z,
								   float &longitude, float &latitude)
{
	latitude = (float)asin(z) * 180.0f / 3.1415927f;
	longitude = (float)atan2(y, x) * 180.0f / 3.1415927f;
}

void globalgrid::latlong_to_threed(float longitude, float latitude,
								   float &x, float &y, float &z)
{
	float lng = longitude*3.1415927f/180.0f;
	float lat = latitude*3.1415927f/180.0f;
	x = (float)(cos(lat) * cos(lng));
	y = (float)(cos(lat) * sin(lng));
	z = (float)sin(lat);
}

void globalgrid::show_coverage(vector<gridcell> &grid)
{
	for (int latitude=0; latitude<(GLOBAL_GRID_DIMENSION/2); latitude++) {
		for (int longitude=0; longitude<GLOBAL_GRID_DIMENSION; longitude++) {
			int y,index = (latitude * GLOBAL_GRID_DIMENSION)+longitude;
			for (y = 1700; y < 3000; y++) {
				if (grid[index].samples[y].size()>0) {
					printf("O");
					break;
				}
			}
			if (y==3000) {
				printf(".");
			}
		}
		printf("\n");
	}
}

void globalgrid::check_grid_coverage(vector<gridcell> &grid)
{
	bool grid_hits[GLOBAL_GRID_DIMENSION*(GLOBAL_GRID_DIMENSION/2)];

	for (int i = 0; i < GLOBAL_GRID_DIMENSION*(GLOBAL_GRID_DIMENSION/2); i++) {
		grid_hits[i]=false;
	}

	for (int latitude=-89; latitude<90; latitude++) {
		for (int longitude=-179; longitude<180; longitude++) {
			int index = get_closest_grid_cell((float)longitude, (float)latitude, grid);
			grid_hits[index] = true;
		}
	}

	int tot=0;
	for (int i = 0; i < GLOBAL_GRID_DIMENSION*(GLOBAL_GRID_DIMENSION/2); i++) {
		if (grid_hits[i]) tot++;
	}

	float coverage = tot*100/(GLOBAL_GRID_DIMENSION*(GLOBAL_GRID_DIMENSION/2));
	printf("Grid centre point coverage with closest cell: %.1f%%\n",coverage);
	if (coverage<50) {

		for (int latitude=0; latitude<(GLOBAL_GRID_DIMENSION/2); latitude++) {
			for (int longitude=0; longitude<GLOBAL_GRID_DIMENSION; longitude++) {
				int index = (latitude * GLOBAL_GRID_DIMENSION)+longitude;
				if (grid_hits[index]) {
					printf("O");
				}
				else {
					printf(".");
				}
			}
			printf("\n");
		}
		printf("\nERROR: Not enough coverage using get_closest_grid_cell\n");
		
	}
}


void globalgrid::check_latlong_conversion(vector<gridcell> &grid)
{
	float longitude2=0,latitude2=0,x=0,y=0,z=0;
	const float tollerance = 2;

	for (int latitude=-89; latitude<90; latitude++) {
		for (int longitude=-179; longitude<180; longitude++) {
			latlong_to_threed(longitude, latitude,
							  x, y, z);
			
			threed_to_latlong(x, y, z,
							  longitude2, latitude2);

			if ((latitude2<-90) || (latitude2>90)) {
				printf("ERROR: Latitude out of range: %.1f\n",latitude2);
			}
			if ((longitude2<-180) || (longitude2>180)) {
				printf("ERROR: Longitude out of range: %.1f\n",longitude2);
			}

			float dlong = (float)longitude - longitude2;
			float dlat = (float)latitude - latitude2;
			if (dlong<0) dlong = -dlong;
			if (dlat<0) dlat = -dlat;

			if (dlong>tollerance) {
				dlong = (float)longitude - -longitude2;
				if (dlong<0) dlong = -dlong;
			}

			if (dlat>tollerance) {
				dlat = (float)latitude - -latitude2;
				if (dlat<0) dlat = -dlat;
			}

			
			if ((dlong>tollerance) || (dlat>tollerance)) {
				if (dlong>tollerance) printf("Longitude error: |%d => %.2f| = %.2f\n",longitude,longitude2,dlong);
				if (dlat>tollerance) printf("Latitude error:  |%d => %.2f| = %.2f\n",latitude,latitude2,dlat);
				printf("ERROR: latitude/longitude conversion is not reversible\n");
				//return;
			}
						
		}
	}
}

// Returns the index of the closest point
int globalgrid:: get_closest_grid_cell(float longitude, float latitude,
									   vector<gridcell> &grid)
{
	float x=0,y=0,z=0,min_dist = 99999;
	int cell_index = 0;

	latlong_to_threed(longitude, latitude, x, y, z);

	for (int i = 0; i < (int)grid.size(); i++) {
		float dx = x - grid[i].x;
		float dy = y - grid[i].y;
		float dz = z - grid[i].z;
		float dist = (dx*dx) + (dy*dy) + (dz*dz);
		if ((i==0) || (dist < min_dist)) {
			cell_index = i;
			min_dist = dist;
		}
	}
	return cell_index;
}

// Sets the centre point of each grid cell
void globalgrid:: points_on_sphere(int no_of_points,
								   vector<gridcell> &grid)
{
    float inc = 3.1415927f * (3.0f - (float)sqrt(5.0));
    float off = 2.0f / no_of_points;
    float x,y,z,r,phi,latitude=0,longitude=0;

    for (int i = 0; i < no_of_points; i++) {
        y = i * off - 1 + (off /2.0f);
        r = (float)sqrt(1 - y * y);
        phi = i * inc;
        x = (float)cos(phi) * r;
        z = (float)sin(phi) * r;

		threed_to_latlong(x, y, z, longitude, latitude);

		grid[i].x = x;
		grid[i].y = y;
		grid[i].z = z;

		grid[i].latitude = latitude;
		grid[i].longitude = longitude;
    }

	check_point_distribution(grid);
}


// Checks the distribution of points on the unit sphere
void globalgrid::check_point_distribution(vector<gridcell> &grid)
{
	float average_separation = 0;
	float min_separation = 99999;
	float max_separation = -999999;

	for (int i = 0; i < (int)grid.size(); i++) {
		float min_dist = 999999;
		bool first = true;
		for (int j = 0; j < (int)grid.size(); j++) {
			if (i!=j) {
				float dx = grid[j].x - grid[i].x;
				float dy = grid[j].y - grid[i].y;
				float dz = grid[j].z - grid[i].z;
				float dist = dx*dx + dy*dy + dz*dz;
				if ((first) || (dist < min_dist)) {
					min_dist = dist;
					first=false;
				}
			}
		}
		average_separation += min_dist;
		if (min_dist < min_separation) min_separation = min_dist;
		if (min_dist > max_separation) max_separation = min_dist;
	}
	average_separation = (float)sqrt(average_separation / (float)grid.size());
	min_separation = (float)sqrt(min_separation);
	max_separation = (float)sqrt(max_separation);

	printf("  Average Grid Cell Separation: %.3f\n",average_separation);
	printf("  Min Grid Cell Separation: %.3f\n", min_separation);
	printf("  Max Grid Cell Separation: %.3f\n", max_separation);
	if (max_separation-min_separation > 0.01) {
		printf("  ERROR: Grid cells not evenly spaced\n");
	}
}

// Generates the grid cells
void globalgrid::generate(vector<gridcell> &grid)
{
	printf("Generating grid\n");

	grid.clear();
	for (int x = 0; x < GLOBAL_GRID_DIMENSION; x++) {
		for (int y = 0; y < GLOBAL_GRID_DIMENSION/2; y++) {
			gridcell cell;
			grid.push_back(cell);
		}
	}

	// set the 3D centre coordinates of each grid cell
	points_on_sphere(GLOBAL_GRID_DIMENSION*(GLOBAL_GRID_DIMENSION/2),grid);

	//check_grid_coverage(grid);
	check_latlong_conversion(grid);
}

void globalgrid::load(vector<tempdata> &data,
					  vector<stationdata> &stations,
					  vector<gridcell> &grid,
					  vector<float> area,
					  vector<long long int> &country_codes,
					  string population_class)
{
	long long int current_station=0;
	float north=0,west=0,altitude=0;
	int latitude_index, longitude_index;
	int grid_cell_index=-1;

	generate(grid);

	printf("Loading anomalies...");
	for (int i = 0; i < (int)data.size(); i++) {
		if (data[i].station!=current_station) {

			bool station_ok=true;
			grid_cell_index = -1;

			if (country_codes.size()>0) {
				station_ok=false;
				for (int c = 0; c < (int)country_codes.size(); c++) {
					if (data[i].country_code == (int)country_codes[c]) {
						station_ok = true;
						break;
					}
				}
			}

			if ((area.size()==4) && (station_ok)) {
			    if (!ghcn::station_within_area(
											   data[i].station,
											   stations,
											   area[0],area[1],
											   area[2],area[3])) {
					station_ok=false;
					current_station = data[i].station;
				}
			}

			if (population_class!="") {
				if (!ghcn::station_has_property(PROPERTY_POPULATION_CLASS,population_class,
												data[i].station,stations)) {
					station_ok=false;
					current_station = data[i].station;
				}
			}

			if (station_ok) {
				if (ghcn::get_station_location(data[i].station,
											   stations,
											   north,west,altitude)) {

					grid_cell_index = get_closest_grid_cell(west,north,grid);
					current_station = data[i].station;
				}
			}
		}
		if (grid_cell_index != -1) {
			grid[grid_cell_index].samples[data[i].year].push_back(data[i]);
		}
	}
	printf("done\n");

	int hits=0;
	for (int i = 0; i < (int)grid.size(); i++) {
		for (int y = 1700; y < 3000; y++) {
			if (grid[i].samples[y].size()>0) {
				hits++;
				break;
			}
		}
	}
	//show_coverage(grid);

	printf("%.1f%% Coverage\n",hits*100.0f/grid.size());
}
